library(dplyr)
# For Model fitting
library(nlme)
library(lme4)
library(purrr)
# For diagnostics
library(performance)
# For adding new columns
library(tibble)
# For highlighting min value in a table
library(reactablefmtr)
# Load data
sys.source("./codes/scripts/code_join_data_full_dataset.R", envir = knitr::knit_global())
# Load functions
sys.source("./codes/functions/functions_mixed_models.R", envir = knitr::knit_global())

Model fixed effects

\[response\sim treatment*fixer\ + initial\ height\ + random(1|specie)\]

Model fitting

# Take response variables' names 
response_vars <- set_names(names(data_complete)[6:(ncol(data_complete))])
model_list <- map(response_vars, ~ mixed_model_1(response = .x, data = data_complete))

Validation plots performance package

map(model_list, check_model)
$rmf


$smf


$lmf


$total_biomass


$above_biomass


$below_biomass


$amax


$gs


$wue


$sla_cm2_g


$perc_c


$perc_n


$ratio_c_n


$d13c


$d15n


$rgr


$agr


$rgr_slope


$Narea_g_m2


$Nmass_mg_g

Model fitting varIdent models for accounting for heterogeneity of varicences

Models that I considered have variance hererogeinity:

varident_response_vars <- response_vars[c("smf","below_biomass","wue","perc_n",
                                          "ratio_c_n","agr","Narea_g_m2","Nmass_mg_g")]

varIdent = None

model_list_base <- map(varident_response_vars, ~ mixed_model_1_lme(response = .x,
                                   data = data_complete))

varIdent = nfixer

model_list_varident <- map(varident_response_vars, ~ mixed_model_2_varident(response = .x,
                                   varident_variable = "nfixer", 
                                   data = data_complete))

varIdent = treatment

model_list_varident_2 <- map(varident_response_vars, ~ mixed_model_2_varident(response = .x,
                                   varident_variable = "treatment", 
                                   data = data_complete))

varIdent = nfixer*treatment

model_list_varident_3 <- map(varident_response_vars, ~ mixed_model_2_varident(response = .x,
                                   varident_variable = "nfixer_treatment", 
                                   data = data_complete))

varIdent = spcode

model_list_varident_4 <- map(varident_response_vars, ~ mixed_model_2_varident(response = .x,
                                   varident_variable = "spcode", 
                                   data = data_complete))

Model comparison AIC values

Validation plots for varIdent models

  • SMF
validation_plot_varident_models("smf", data = data_complete, 
                                varident_variable = "spcode")

  • below_biomass
validation_plot_varident_models("below_biomass", data = data_complete, 
                                varident_variable = "nfixer_treatment")

  • wue
validation_plot_varident_models("wue", data = data_complete, 
                                varident_variable = "nfixer_treatment")

  • perc_n
validation_plot_varident_models("perc_n", data = data_complete, 
                                varident_variable = "spcode")

  • ratio_c_n
validation_plot_varident_models("ratio_c_n", data = data_complete, 
                                varident_variable = "spcode")

  • agr
validation_plot_varident_models("agr", data = data_complete, 
                                varident_variable = "spcode")

  • Narea_g_m2
validation_plot_varident_models("Narea_g_m2", data = data_complete, 
                                varident_variable = "spcode")

  • Nmass_mg_g
validation_plot_varident_models("Nmass_mg_g", data = data_complete, 
                                varident_variable = "spcode")

Models selected